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interaction  fields  of  a shock  wave  with  a boundary  layer,  both 
for  the  laminar  and  for  the  turbulent  cases  with  and  without 
separation.  The  latest  algorithm  is  developed  specifically  for 
computational  solutions  at  coarse  mesh  attainable  with  a modest 
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Princeton  University.  It  is  shown  how  these  solutions  with  com- 
parable accuracy  can  be  achieved  with  an  order  of  magnitude 
reduction  in  computer  time  and  capability  when  the  new  algorithm 
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FINAL  REPORT  TO  OFFICE  OF  NAVAL  RESEARCH 
Contract  No.  N 00014-75-C-0376 
"The  Application  of  Difference  Methods  in  Fluid  Mechanics" 


The  primary  objective  of  the  proposed  research  has  been 
the  development  of  difference  methods  for  the  solution  of  gas 
dynamic  equations  as  a tool  of  fluid  dynamics  research  and 
developmental  work.  Fundamental  understanding  is  emphasized 
rather  than  the  generation  of  specific  solutions  of  various 
practical  flow  problems;  although  it  has  been  conducted  around 
two  specific  problems. 

The  first  problem  is  the  flow  field  of  interaction  of  a 
shock  wave  with  a planar  boundary  layer.  The  second  problem  is 
the  supersonic  flow  over  a yawed  cone.  The  first  problem  is 
distinctly  in  two  space  dimensions  while  the  second  is  dominated 
by  important  three  dimensional  effects. 

I I ' 


The  requirement  of  adequate  resolution  for  the  computational 
solution  of  3D  flow  problems  makes  it  necessary  to  consider  the 
use  of  the  very  large  and  very  high  speed  computers  like  the 
ASC  or  the  STAR  machines.  The  ASC  computer  was  operational  and 
available  in  1977.  In  earlier  years  a certain  amount  of 
analysis  was  done  and  a preliminary  computation  was  performed 
on  the  IBM  360-91  computer  to  learn  about  the  various  subtlties 
of  the  computational  solution  of  the  truely  3-D  flows.  We 
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began  constructing  the  3-D  program  for  the  computational  solu- 
tion of  the  yawed  cone  problem  working  on  the  ASC  machine  in 
1977.  We  have  low  priority  for  our  access  to  the  ASC  machine 
and  have  progressed  by  early  1978  only  to  the  successful  test 
of  the  unyawed  configuration.  It  was  decided  in  1978  as  was 
reported  to  the  ONR  office  in  Washington  D.C.  to  devote  our 
effort  to  the  first  problem  in  two  space  dimensions  in  the  re- 
maining period.  Our  schemes  conceived  earlier  for  generating 
the  yawed  cone  solution  from  the  axisymmetric  situation  was  not 
tested  on  the  ASC  machine.  A relatively  complete  picture  is  ob- 
tained for  the  first  problem  as  a result. 

From  the  preliminary  studies  of  3-D  problems,  the  following 
"tentative  conclusions"  were  obtained. 

1.  The  stability  of  our  computational  algorithm  is  only  slightly 
affected  in  the  interior  (of  the  field  of  computation) . The  ranges 
of  Mesh  Reynolds  number  for  stable  computations  remains  essenr 
tially  unchanged,  although  the  maximum  amplitudes  of  oscillatory 
perturbations  that  permit  stable  computations  can  be  significantly 
reduced  at  certain  ranges  of  values  of  Mesh  Reynolds  numbers. 

The  reduction  of  such  maximal  allowable  disturbance  amplitudes 
is  really  not  harmful,  since  the  lower  limits  of  the  amplitudes 
are  too  large  to  be  tolerated  in  a reasonable  solution  anyway. 

2.  Proper  treatment  of  the  initial  boundary  conditions  is  con- 
siderably more  critical  in  3-D  then  in  2-D  problems  especially 
where  shockwaves  intersect  the  computational  boundary.  This 
aspect  can  be  very  important  in  the  computational  solution  of 
such  a complex  problem  as  the  supersonic  flow  over  a yawed  cone. 
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We  have  had  no  opportunity  to  try  out  in  the  full  3-D  field  any 
schemes  of  boundary  treatment  analyzed  in  model  study. 

It  was  decided  to  suspend  the  investigation  of  the  yawed 
core  problem  in  the  spring  of  1978,  although  certain  analytic 
aspects  of  the  yawed  cone  problem  continued  until  fall.  No 
concrete  results  can  be  reported.  Partial  documentation  of  the 
above  results  are  contained  in  the  two  Master  Theses  by  Sung 
and  Specht  and  the  incomplete  manuscript  of  Mr.  Grasso. 

The  treatment  of  the  shock  wave  boundary  layer  interaction 
problem  was  divided  into  three  stages: 

1.  Development  of  algrithms  and  their  preliminary  analysis  and 
tests.  The  test  computations  include  those  for  shock  waves  in 
one  space  dimension  and  then  in  two  space  dimensions. 

2.  Computational  solution  of  the  interaction  field  resulting  . 
from  the  incidence  of  a weak  shock  wave  on  a laminar  boundary 
layer  with  and  without  separation. 

3.  Computational  solution  of  the  interaction  field  resulting 
from  the  incidence  of  a moderate  shock  on  a turbulent  boundary 
layer  with  separation. 

It  is  generally  presumed  that  by  computing  at  sufficiently 
small  meshes,  the  computed  results  would  be"sufficiently  accu- 
rate" although  the  errors  may  not  be  commeasurate  with  the  local 
truncation  errors.  Such  is  true  and  proven  for  those  linear, 
smooth,  well  posed  problems.  For  those  practical  problems,  which 
are  nonlinear,  not  so  smooth,  and  often  poorly  posed,  this  pre- 
sumption is  at  best  a plausible  conjecture.  This  was  discussed 
and  exemplified  quite  extensively  in  the  VKI  lecture  notes  (1974) 
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by  the  author  [4], 

Our  originial  approach  to  the  solution  of  the  shock  wave 
laminar  boundary  layer  interaction  problem  was  based  on  such 
conventional  wisdom  of  successive  mesh  refinements;  an  approach 
which  was  apparently  rather  successful  in  our  earlier  inves- 
tigation of  similarly  complex  problems  like  near  wake  and 
leading  edge  flows.  We  did  obtain  very  reasonable  solutions 
(Ph.D.  Thesis (Messina) ) , in  much  better  agreement  with  experi- 
mental data  than  other  available  computational  solutions. 

The  behavior  of  the  successive  iterates  of  Messina's  solu- 
tion raises  however,  some  serious  questions  as  was  suggested 
by  our  linearized  analysis  of  the  phase  errors  in  the  solution. 
The  oscillatory  error  components  in  the  computed  solution  tends 
to  be  minimal  when  the  phase  of  these  oscillations  change 
sign.  Accordingly  the  error  in  the  computational  solution  may 
reach  a local  minimum  at  some  finite  value  of  mesh  Reynolds 
number.  Mesh  refinement  beyond  the  critical  value  can  cause 
deterioration  of  the  computed  solution. 

The  above  preliminary  result  is  amplified  and  confirmed  by 
the  detailed  analysis  of  the  system  of  nonlinear  difference 
equations  of  one  dimensional  model  equations.  This  analytic 
study  (Ph.D.  Thesis (Shubin) ) , confirmed  the  asymptotic  nature 
of  the  convergence  of  the  successive  iterants  with  mesh  re- 
finements i.e.  reducing  the  mesh  Reynolds  number.  The  con- 
vergence is  likely  to  be  non-uniform,  especially  when  the  dif- 
ference formulation  includes  those  extraneous  boundary  condi- 
tions commonly  adopted  in  the  solution  of  practical  problems. 
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The  errors  in  the  computed  solution  may  actually  increase  to 
fairly  large  values  with  further  mesh  refinements  beyond  the 
optimal  limit.  Whether  the  sequence  of  the  iterants  will 
eventually  converge  to  the  correct  solution  as  Re^x  0 will 
depend  on  how  well  the  difference  problem  is  posed,  whether  the 
formulation ?is  consistent  and  stable,  and  may  be  on  other 
conditions.  But  it  is  apparent  that  when  most  common  differ- 
ence schemes  are  applied  to  complex  practical  problems,  such  a 
convergence  region  of  sufficiently  small  Re^x  even  if  it  exists, 
is  well  beyond  the  reach  of  computers  in  the  forseeable  future. 

The  analysis  indicates  further,  that  the  magnitudes  of  the 
minimal  error  of  the  computed  solutions  at  the  critical  mesh 
Reynolds  number  can  be  "comparable"  to  that  of  the  convergent 
sequence  at  much  smaller  mesh  Reynolds  numbers. Thus,  such 
computed  solutions  near  the  minimal  errors  at  fairly  large  Re^x 
may  be  as  useful  as  the  convergent  approximations  at  much 
smaller  Re^x.  The  situation  is  quite  similar  to  the  approximate 
evaluation  of  transcendental  functions  as  the  solutions  of  some 
ordinary  differential  equations.  Such  asymptotic  approximations 
are  then  much  more  useful  in  practice  than  those  convergent  ap- 
proximations. 

When  the  class  of  simple  difference  algorithms  which  we  anal- 
yzed (Shubin) , is  applied  to  complex  practical  problems  in  multi- 
space dimensions,  there  is  an  outstanding  difficulty.  The  mesh 
Reynolds  number  generally  varies  over  a much  wider  range  of  values 
in  a complex  flow  problem  than  the  plateau  of  errors  around  the 
critical  mesh  Reynolds  number.  Moreover,  the  stability  of  the  com- 
plete difference  formulation  will  depend 
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on  the  discretization  parameter  y in  the  class  of  algorithm 
analyzed  and  there  are  indications  of  computational  instability 
if  y is  too  large  negatively.  In  order  to  take  advantage  of 
such  asymptotic  approximations  of  coarse  mesh  computational 
solution,  we  have  to  develope  some  algorithm  that  promises 
stable  computation  and  small  and  very  flat  variations  of 
errors  at  large  mesh  Reynolds  numbers.  This  is  accomplished 
by  a simple  modification  of  our  two  step  algorithms  . 

Both  the  predictor  and  the  correctors  of  the  two  step 
scheme  of  Cheng-Alien  (Scheme  1)  are  defined  over  a single 
mesh  (or  grid)  system.  We  now  consider  the  predictor  step  de- 
fined on  a sub-grid  at  mid-mesh  points  and  let  the  function 
values  at  mid-mesh  be  the  arithematic  averages  of  the  nearest 
neighbors.  This  new  scheme  (2)  is  obtained  in  a manner  simi- 
lar to  Richtmyer's  modification  of  the  Lax-Wendroff  scheme. 
Linearized  analysis  provides  very  encouraging  signs  as  a scheme 
particularly  suitable  for  coarse  mesh  computation.  The  scheme 
was  further  tested  computationally  with  nonlinear  Burgers' 
equation  and  then  the  Navier-Stokes  equations  in  two  space 
dimension  for  the  simple  flow  field  resulting  from  the  deflec- 
tion of  a uniform  supersonic  stream  by  a straight  oblique  shock. 
The  following  encouraging  results  are  obtained: 

1.  The  maximal  errors  of  the  computed  solution  with  scheme  (2) 
is  significantly  below  those  with  scheme  (1)  at  large  mesh 
Reynolds  numbers  (ReAx  £ 3) , although  the  reverse  is  true  for 


Reix  * 3. 


2.  Over  the  entire  range  of  Re^x  scheme  2 appears  to  provide 
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solutions  well  within  the  accuracy  reauiremsnts  (<  10%)  in 
engineering  applications,  and  comparable  to  the  accuracy  of  the 
data. 

3.  The  oscillations  (or  rather  overshoots)  in  results  computed 
with  scheme  2,  if  they  should  occur, are  confined  within  a mesh 
or  two  away  from  the  shock  like  discontinuity. 

4.  Scheme  2 converges  to  the  "steady  state"  much  faster  than 
scheme  1,  i.e.  fewer  number  of  steps  are  required  to  satisfy  a 
given  steady  state  criterion  from  a given  initial  state. 

The  same  program  for  the  Navier-Stokes  equation  is  then 
used  to  solve  the  interaction  field  resulting  from  the  incidence 
of  a weak  shock  wave  on  a laminar  boundary  layer.  In  recogni- 
tion of  the  nontrival  pressure  gradient  along  the  wall  and  of 
the  coarseness  of  our  mesh,  parabolic  mean  velocity  profiles 
is  used  in  the  treatment  of  the  wall  cells.  The  cases  to  be 
calculated  correspond  to  the  experimental  data  of  Hakkinen  et. 
al.  This  unseparated  case  has  been  computationally  solved  by 
many  authors  including  our  solution  by  Messina  with  scheme  1, 
which  represents  the  most  "accurate"  computational  solution  as 
judged  by  comparison  with  experimental  data.  Our  interest  is 
to  demonstrate  that  with  scheme  2,  solutions  as  "accurate"  as 
the  above  can  be  obtained  at  much  coarser  mesh,  and  with  much 
less  demand  on  computing  time  and  computer  capability.  As  is 
illustrated  in  the  following  table  (see  next  page) , our  ob- 
jectives are  largely  met  with  1/8  of  the  number  of  grid  points 
in  the  solution  with  scheme  2,  the  computer  time  required  to 
reach  steady  state  from  the  same  initial  data  is  reduced  by  a 


factor  1/30.  The  200  K-bytes  of  the  computer  storage  required 
is  largely  to  accomodate  the  program  (175  K-bytes) , which,  if 
streamlined,  can  be  substantially  reduced.  Thus  a problem  of 
comparable  complexity  can  be  handled  with  computers  with  a 
fraction  of  the  capabilities  of  the  IBM  360-91.  Using  IBM  360- 
91  it  is  possible  to  compute  somewhat  more  complex  problems. 

The  separated  case  of  Hakkinen's  data  was  not  successfully 
computed  with  scheme  1,  possibly  because  of  poor  stability  and/ 
or  slow  convergence  toward  steady  state.  This  case  is  now 
solved  with  scheme  2 at  half  the. mesh  size;  and  the  steady  state 
criterion  was  satisfied  after  1200  cycles  (vs.  600  for  un- 
separated case)  and  15  minutes  (v.s.  2.4  minutes)  of  computation 
with  IBM  360-91.  The  results  compare  as  favorable  with  ex- 
perimental data  as  are  the  unseparated  case.  The  separation 
bubble  is  very  substantial,  which  makes  it  necessary  to  extend 
the  field  of  computation  to  766^  compared  with  586 ^ for  the 
unseparated  case. 

The  successful  completion  of  the  separated  case  of  shock- 
wave  laminar  boundary  layer  interaction  encourages  us.  to 
obtain  computational  solution  of  similar  interaction  problem 
with  turbulent  boundary  layer.  There  remains  two  outstanding 
problems.  The  first  one  is  physical  that  a model  for  the 
turbulence  stress  is  needed.  The  second  one  is  computational 
since  it  is  as  yet  uncertain  whether  scheme  2 at  coarse  mesh 
can  successfully  handle  the  "intersecting  shockwaves"  i.e. 
the  incident  shock  and  the  induced  shock  wave  sprang  from  the 
"sudden"  response  of  the  turbulent  boundary  layer  to  downstream 
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pressure  rise. 

To  focus  our  attention  on  the  computational  aspects  of  the 
problem,  we  decided  to  adopt  an  eddy  diffusivity  model  of  tur- 
bulent stress  so  that  the  governing  differential  equations 
system  will  still  be  the  Navier-Stokes . The  eddy  diffusivity 
e is  chosen  algebraically  so  as  to  represent  the  stress 
variation  across  the  turbulent  boundary  layer  properly.  Thus 
in  the  outer  region  e is  specified  according  to  clauser,  and 
in  the  inner  wall  region  to  the  model  of  Townsend  and  Mellor. 
Since  separation  is  present  in  the  case  to  be  calculated,  the 
wall  model  was  extended  further  to  include  separation  and 
separated  layers. 

Computational  solutions  were  then  attempted  for  the  following 

cases  in  a Mach  3 supersonic  stream  with  shocks,  giving  the 

flow  deflection  angles  0 = 7.93°,  9.87°  and  10.83° , incident  on 

a turbulent  boundary  layer  with  Re.,  - 1.5  x 10s.  The  cases 

0 1 

correspond  to  the  experimental  data  of  Law  and  have  also  been 
solved  computationally  by  Shank  et.al.  Shank  et.al.  employed 
a slightly  different  eddy  diffusivity  model  and  introduced  a 
relaxation  parameter  X which  was  adjusted  to  achieve  fit  with 
known  experimental  results  of  Law.  The  present  results  for 
these  cases  agree  with  Law's  experimental  data  as  well  but 
without  the  adjustable  relaxation  parameter.  All  cases  involve 
some  region  of  separated  flow. 

One  more  case  with  0 = 12.75  at  Re.,  = 2.5  x 10s  was  com- 

o l 

puted.  It  corresponds  to  the  experimental  data  of  Bogdonoff 
and  Kepler.  This  case  contains  the  largest  separated  region. 
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It  has  been  computed  by  Wilcox  based  on  Saffman's  differential 
model  of  turbulence.  The  present  results  with  eddy  diffusivity 
compares  more  favorably  with  the  experimental  data.  There  is 
clearly  room  for  improvement  especially  quantitatively;  but 
the  computed  results  appear  quite  satisfactory.  It  is  not 
possible  to  identify  how  much  of  the  existing  disparity  between 
the  experimental  data  and  the  computed  results  is  due  to  errors 
in  the  computational  procedure  and/or  those  in  the  turbulence 
model.  We  are  particularly  encouraged  by  the  fact  that  our 
coarse  mesh  scheme  can  compute  successfully  the  configuration 
of  intersecting  shocks  as  well  as  the  abrupt  generation  of  a 
shock  wave  from  the  sudden  changes  in  the  turbulent  boundary 
layer.  Both  the  two  shock  and  the  three  shock  interaction  con- 
figuration have  been  successfully  computed.  The  use  of  coarse 
mesh  and  the  relatively  fast  convergence  with  the  present 
scheme  2 permits  significant  savings  in  computer  time.  This  is 
illustrated  in  the  following  Table. 

In  conclusion,  we  have  developed  a new  difference  scheme 
(scheme  2)  particularly  suited  for  the  computational  solution 
of  complex  flow  problems  with  rather  coarse  mesh.  The  various 
examples  of  the  computational  solution  of  the  interaction  field 
resulting  from  the  incidence  of  a shock  wave  on  a boundary  layer 
with  or  without  separation,  laminary  or  turbulent,  are  probably 
as  complex  as  many  problems  of  practical  interest.  The  new 
scheme  is  simple  and  convenient  to  use.  It  possesses  good 
stability  property  and  appears  to  converge  toward  the  "steady 
state"  considerably  faster  than  our  earlier  scheme  1.  It  is 
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particularly  suitable  for  coarse  mesh  computations.  All  these 

I* 

factors  contribute  to  significant  savings  in  computer  time. 

Moreover,  it  opens  up  the  possibility  of  the  computational.' 
solution  of  rather  complex  flow  problems  with  modest  computer 
requirements.  We  have  thus  reached  a definite  stage  of  developing 

i 1 : 

difference  methods  for  the  solution  of  gas  dynamic  equations 
as  a tool  of  fluid  dynamics  research  and  development.  Details 
of  the  development  and  test  of  the  new  algorithms  are  largely 
contained  in  the  Ph.D.  thesis  of  L.  Oey.  Technical  reports 
are  being  prepared  for  submission  for  Journal  publication. 

They  will  be  distributed  according  to  the  ONR  distribution 
list  when  published. 
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Documentation  Generated 
Journal  Publications 


During  the  contractural  period  1972  - 1979,  the  following 
technical  publications  have  been  generated,  or  are  under 
preparation. 


1.  "Finite  Difference  Treatment  of  Strong  Shock  over  a sharp 
Leading  Edge  with  Navier-Stokes  Equations".  (with  J.H. 

Chen) , Proceedings  of  3rd  International  Conference  on 
Numerical  Methods  in  Fluid  Mechanics.,  Paris  (1972), 
Published  in  Lecture  Notes  in  Physics  Series  No.  18  (1973) , 
Springer  Verlag. 

2.  "Friction  and  Heat  Transfer  Laws  in  Slip  Flows",  Proceedings 
in  10th  International  Symposium  on  Space  Technology  and 
Science,  Tokyo,  Japan  1973. 

3.  "Slip  Friction  and  Heat  Transfer  Laws  in  a Merged  Region" 
(with  J.H.  Chen),  The  Physics  of  Fluid,  Vol,  17,  No.  9 
(1974)  . 

4.  "A  Critical  Review  of  Numerical  Solutions  of  Navier-Stokes 
Equations",  Lecture  series  delivered  at  Von  Nauman 
Institute  of  Aerodynamics,  Belgium,  1974.  Published  in 
"Progress  in  Numerical  Fluid  Dynamics",  Edited  by  H.J. 

Wirz,  Lecture  note  in  Physics  series  No.  41  (1975)  Springer 
Verlag. 

5.  "Computational  Accuracy  and  Mesh  Reynolds  Number",  (with 

.G.  Shubin) , Journal  of  Computational  Physics  2£  (1978) . 

6.  "One  Dimensional  Gas  Dynamics  Model  Study  of  Computational 
Accuracy"  ( .G.  Shubin) , accepted  to  appear  in  Journal  of 
Computational  Physics. 

7.  "Errors  in  Finite  Difference  Solutions  of  Navier-Stokes 
Equations",  Proceedings  of  6th  International  Conference  on 
Numerical  Methods  in  Fluid  Mechanics.  Tbilis  U.S.S.R. 
(1978)  to  be  published  in  Lecture  Notes  in  Physics  Series. 
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"Coarse  Mesh  Computational  Solution  of  Navier-Stokes 
Equations",  (with  L.  Oey)  in  preparation,  to  be  submitted 
for  Journal  publication. 


...  _ 


L3 


"Interaction  of  Shock  Wave  with  Laminar  Boundary  Layer" 
to  be  prepared  with  (N.  Messina  and  L.  Oey) . 

"Computational  Solution  of  the  Interaction  of  a Shock 
Wave  with  a Turbulent  Boundary  Layer" , to  be  prepared 
with  (L.  Oey). 
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University  Reports 
Ph.D.  Thesis 


1.  J.H.  Chen,  "Finite  Difference  Methods  and  the  Leading  Edge 
Problem",  1971. 

2.  N.A.  Messina,  "A  Computational  Investigation  of  Shock  Waves, 
Laminar  Boundary  Layers  and  their  Mutual  Interaction",  1977. 

3.  G.R.  Shubin,  "One  Dimensional  Gas  Dyanmic  Modeling  and 
Computational  Accuracy",  1977. 

4.  L.Y.  Oey,  "Large  Mesh  Reynolds  Number  Computations  - with 
Applications  to  Shock  Boundary  Layer  Interaction  Problems", 
1979. 

5.  K.  Meintjes,  in  preparation,  expected  1979 
(.tentative  title  - "Predictor-Corrector  Methods  for  Time 
Dependent  Compressible  Flows"). 

6.  F.  Grasso,  "Computational  Solution  of  Supersonic  Flow  over 
a Yawed  Core".  (incomplete). 


